Introdução a Mapas no R

Nicholas A. C. Marino

Elementos da Aula (baseada na aula da Dra. Andrea Sanchez Tapia/2016) + atualizações minhas

1. O que está aula é e o que ela não é!;
2. Utilizando o conteúdo do R para gerar mapas simples;
3. Criando mapas simples a partir de arquivos raster;
4. Usando informações da interação com a internet para criar mapas simples.

O que está aula é e o que ela não é!

O R possui uma ampla capacidade e versatilidade para a execução de praticamente qualquer tarefa - inclusive para a criação e edição de mapas. É possível usar o R e toda a sua funcionalidade de forma muito similar àquela de programa como o ArcGIS. No entanto, este não será o nosso objetivo aqui, mas sim apresentar a funcionalidade mais básica possível do R para a criação de mapas (principalmente aqueles que sempre queremos adicionar em uma dissertação, tese ou manuscrito).

Existe uma diversidade enorme de material e sites na internet com material sobre a confecção de mapas no R. Sugiro dar uma olhada de forma mais detalhada no material gerado pela Dra. Andrea Sanchez Tapia para este mesmo curso em 2016 (99_cursos_anteriores/09_Rmaps) e/ou nos sites listados abaixo:

• R Special Interest Group on using Geographical data and Mapping: https://stat.ethz.ch/mailman/listinfo/r-sig-geo
  
• TaskView: Analysis of Spatial Data: https://cran.r-project.org/web/views/Spatial.html
  
• Lovelace, R., & Cheshire, J. (2014). Introduction to visualising spatial data in R. National Centre for Research Methods Working Papers, 14(03). Retrieved from https://github.com/Robinlovelace/Creating-maps-in-R
• Kahle & Wickham 2012 - ggmaps: https://dl.dropboxusercontent.com/u/24648660/ggmap%20useR%202012.pdf
  
• Kim Gilbert http://www.zoology.ubc.ca/~kgilbert/mysite/Miscellaneous_files/R_MakingMaps.pdf
  
• Kevin Johnson http://www.kevjohnson.org/making-maps-in-r/
  
• Santiago Begueria. Mapping with ggplot2: hexbin maps. http://santiago.begueria.es/2016/04/mapping-with-ggplot2-hexbin-maps/
  
• http://geog.uoregon.edu/GeogR/topics/maps01.html
  
• A função para acrescentar escala em ggplot2: http://editerna.free.fr/wp/?p=76

Alguns pacotes úteis para a confecção de mapas no R são:

• RgoogleMaps: baixar mapas do Google;
  
• raster: pacote para a criação e manipulação de arquivos raster;
  
• maptools: ler ESRI shapefiles;
  
• rgdal;
  
• maps: mapas simples, eixos, escalas e cidades;
  
• mapdata: base de dados WorldHires e rios;
  
• rworldmap: outra base de dados de mapas do mundo;
• ggplot2;
  
• ggmap: mapas customizados, open street maps e outras funcionalidades de mapas para o ggplot2;
  
• ggalt: outras funcionalidades para a customização de mapas e escalas;
  
• rasterVis: para facilitar a integração de arquivos raster com o ggplot2;
  
• …

Utilizando o conteúdo do R para gerar mapas simples

O pacote maps pode plotar mapas simples através de dados existentes no próprio pacote:

# carregando pacotes
library(maps) # plotagem genérica de mapas
library(mapdata) # base de dados worldHires
# plotando o mapa do Brasil
map(database = "worldHires", regions = "Brazil")

Você pode adicionar a projeção dos outros países no espaço disponível:

map(database = "worldHires", regions = "Brazil")
map(add = TRUE)

Adicionando os eixos de longitude e latitude:

map(database = "worldHires", regions = "Brazil")
map(add = TRUE)
map.axes()

Adicionando escala ao mapa e uma linha para representar o Equador:

map(database = "worldHires", regions = "Brazil")
map(add = TRUE)
map.axes()
map.scale(ratio = FALSE, cex = 0.7)
abline(h = 0, lty = 2)

Você também pode plotar os pontos correspondentes às cidades ao mapa usando a função:

map(database = "worldHires", regions = "Brazil")
map(add = TRUE)
map.axes()
map.scale(ratio = FALSE, cex = 0.7)
abline(h = 0, lty = 2)
map.cities(country = "Brazil", minpop = 1000000, pch = 19, cex = 1.2)

---

Exercício 1

Modifique as opções apresentadas anteriormente para explorar as opções gráficas disponíveis nas funções utilizadas.

---

Uma outra opção para criar mapas com base nos dados já existentes nos pacotes do R é através do ggplot2, através da função borders. A função borders (do ggplot2) funciona de forma similar àquela do pacote maps.

# carregando o ggplot2
suppressPackageStartupMessages(library(ggplot2))

---

Exercício 2

Utilize a função borders do pacote ggplot2 para criar um mapa editado do Brasil.

---

Criando mapas simples a partir de rasters e shapefiles

Como vimos alguns dos pacotes disponíveis no R já oferecem mapas embutidos em suas funções. Todavia, muitas destas funções e seus mapas não satisfazem totalmente algumas das necessidades que alguns podem ter - e outros não satisfazem nem um pouco. Na realidade, na maior parte das vezes, é mais útil termos um conjunto de shapefiles ou de rasters a partir dos quais vamos conseguir construir os mapas que de fatos queremos.

Você pode encontrar arquivos para a confecção de mapas em:

• Shapefiles: pontos, linhas e polígonos (e.g., cidades, rios, estados…)
  • do Brasil: ftp://geoftp.ibge.gov.br/mapas_interativos/
  • do mundo: DIVA_GIS project http://www.diva-gis.org/Data e por país: http://www.diva-gis.org/gdata
  • Global Administrative Areas http://gadm.org (você pode baixar estes dados diretamente do R através da função raster::getData)
  • Google!
• Rasters: superfícies contínuas, reticuladas (grids), com diferentes tamanhos de pixel.
  • modelos digitais de elevação: http://eros.usgs.gov/
  • clima: Worldclim http://www.worldclim.org/ (você pode baixar estes dados diretamente através da função raster::getData)

A função raster::getData é particularmente útil para baixarmos os rasters e shapefiles diretamente de base de dados online. Não vamos fazer o download de um conjunto de rasters aqui pois necessitamos que a internet esteja totalmente funcional; todavia, podemos observar um raster que já foi baixado a partir do “worldclim” abaixo:

# carregando o pacote raster
suppressPackageStartupMessages(library(raster))
# carregando o arquivo de raster do worldclim
bio1 <- raster("layers/bio1.bil")
# plotando o raster
plot(bio1, xlab = "Longitude", ylab = "Latitude", main = "Temperature Média do Ar (x 10 ºC)")

Você pode descobrir as coordenadas de uma localidade usando a função ggmap::geocode, como se estivesse fazendo uma busca nos mapas do Google.

# carregando pacote
suppressPackageStartupMessages(library(ggmap))
# quais as coordenadas da ilha do fundão?
geocode(location = "Ilha do Fundão")

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Ilha%20do%20Fund%C3%A3o&sensor=false

##        lon       lat
## 1 -43.2259 -22.86269

Você pode extrair os dados das coordenadas de um objeto do tipo raster usando a função raster::extract, que pede um objeto raster x e uma matriz ou data.frame como um objeto y:

# temperatura média anual do ar na ilha do fundão (x 10ºC)
extract(x = bio1, y = cbind(-43.2259, -22.86269))

##     
## 224

Para facilitar a manipulação de arquivos com formato raster, você pode usar o pacote rasterVis e a função gplot para plotar o arquivo raster utilizando as funções do ggplot.

# carregando pacote
suppressPackageStartupMessages(library(rasterVis))
# plotando figura
gplot(bio1) +
  geom_tile(aes(fill = value)) +
  scale_fill_gradient(low = 'blue2', high = 'firebrick1') +
  coord_equal()

---

Exercício 3

1. Carregue o arquivo raster bio12.bil, que corresponde aos dados da precipitação anual média.
   
2. Extraia as coordenadas de uma localidade de sua escolha.
   
3. Determine a precipitação anual média para a localidade de sua escolha.
4. Crie e customize um gráfico utilizando este arquivo raster.

---

Além de arquivos raster, podemos carregar shapefiles no R através da função readShapeSpatial do pacote maptools (ou variantes: readShapeLines, readShapePoints, readShapePoly). Por exemplo:

# carregando pacote
library(maptools)
# carregando o shapefile
brasil <- readShapeSpatial("layers/BRA_adm/BRA_adm3.shp")

O shapefile é um objeto que também pode ser manipulado e plotado individualmente!

# somente as dez primeiras linhas
brasil@data[1:10,]

##   ID_0 ISO NAME_0 ID_1 NAME_1 ID_2          NAME_2 ID_3          NAME_3
## 0   33 BRA Brazil    1   Acre    1      Acrelândia    1      Acrelândia
## 1   33 BRA Brazil    1   Acre    2    Assis Brazil    2    Assis Brazil
## 2   33 BRA Brazil    1   Acre    3       Brasiléia    3       Brazileia
## 3   33 BRA Brazil    1   Acre    4          Bujari    4          Bujari
## 4   33 BRA Brazil    1   Acre    5        Capixaba    5        Capixaba
## 5   33 BRA Brazil    1   Acre    6 Cruzeiro do Sul    6 Cruzeiro do Sul
## 6   33 BRA Brazil    1   Acre    7  Epitaciolândia    7  Epitaciolândia
## 7   33 BRA Brazil    1   Acre    8           Feijó    8           Feijó
## 8   33 BRA Brazil    1   Acre    9          Jordão    9          Jordão
## 9   33 BRA Brazil    1   Acre   10     Mâncio Lima   10     Mâncio Lima
##     TYPE_3 ENGTYPE_3 NL_NAME_3 VARNAME_3
## 0 Distrito  District      <NA>      <NA>
## 1 Distrito  District      <NA>      <NA>
## 2 Distrito  District      <NA>      <NA>
## 3 Distrito  District      <NA>      <NA>
## 4 Distrito  District      <NA>      <NA>
## 5 Distrito  District      <NA>      <NA>
## 6 Distrito  District      <NA>      <NA>
## 7 Distrito  District      <NA>      <NA>
## 8 Distrito  District      <NA>      <NA>
## 9 Distrito  District      <NA>      <NA>

# cortando o shapefile para o estado
estado_rio <- subset(brasil, brasil@data$NAME_1 == "Rio de Janeiro")
# observando o shapefile
estado_rio@data[1:10,]

##      ID_0 ISO NAME_0 ID_1         NAME_1 ID_2         NAME_2 ID_3
## 6683   33 BRA Brazil   19 Rio de Janeiro 3564 Angra dos Reis 6684
## 6684   33 BRA Brazil   19 Rio de Janeiro 3564 Angra dos Reis 6685
## 6685   33 BRA Brazil   19 Rio de Janeiro 3564 Angra dos Reis 6686
## 6686   33 BRA Brazil   19 Rio de Janeiro 3564 Angra dos Reis 6687
## 6687   33 BRA Brazil   19 Rio de Janeiro 3564 Angra dos Reis 6688
## 6688   33 BRA Brazil   19 Rio de Janeiro 3564 Angra dos Reis 6689
## 6689   33 BRA Brazil   19 Rio de Janeiro 3565        Aperibé 6690
## 6690   33 BRA Brazil   19 Rio de Janeiro 3566       Araruama 6691
## 6691   33 BRA Brazil   19 Rio de Janeiro 3566       Araruama 6692
## 6692   33 BRA Brazil   19 Rio de Janeiro 3566       Araruama 6693
##                    NAME_3   TYPE_3 ENGTYPE_3 NL_NAME_3 VARNAME_3
## 6683               Abraao Distrito  District      <NA>      <NA>
## 6684       Angra dos Reis Distrito  District      <NA>      <NA>
## 6685           Cunhambebe Distrito  District      <NA>      <NA>
## 6686           Jacuecanga Distrito  District      <NA>      <NA>
## 6687            Mambucaba Distrito  District      <NA>      <NA>
## 6688   Praia de Aracatiba Distrito  District      <NA>      <NA>
## 6689              Aperibé Distrito  District      <NA>      <NA>
## 6690             Araruama Distrito  District      <NA>      <NA>
## 6691         Morro Grande Distrito  District      <NA>      <NA>
## 6692 São Vicente de Paula Distrito  District      <NA>      <NA>

# cortando o shapefile para a cidade
cidade_rio <- subset(brasil, brasil@data$NAME_2 == "Rio de Janeiro")
# observando o shapefile
cidade_rio@data

##      ID_0 ISO NAME_0 ID_1         NAME_1 ID_2         NAME_2 ID_3
## 6929   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6930
## 6930   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6931
## 6931   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6932
## 6932   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6933
## 6933   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6934
## 6934   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6935
## 6935   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6936
## 6936   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6937
## 6937   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6938
## 6938   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6939
## 6939   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6940
## 6940   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6941
## 6941   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6942
## 6942   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6943
## 6943   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6944
## 6944   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6945
## 6945   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6946
## 6946   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6947
## 6947   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6948
## 6948   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6949
## 6949   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6950
## 6950   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6951
## 6951   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6952
## 6952   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6953
## 6953   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6954
## 6954   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6955
## 6955   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6956
## 6956   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6957
## 6957   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6958
## 6958   33 BRA Brazil   19 Rio de Janeiro 3630 Rio de Janeiro 6959
##                  NAME_3   TYPE_3 ENGTYPE_3 NL_NAME_3 VARNAME_3
## 6929           Anchieta Distrito  District      <NA>      <NA>
## 6930              Bangu Distrito  District      <NA>      <NA>
## 6931    Barra da Tijuca Distrito  District      <NA>      <NA>
## 6932           Botafogo Distrito  District      <NA>      <NA>
## 6933       Campo Grande Distrito  District      <NA>      <NA>
## 6934             Centro Distrito  District      <NA>      <NA>
## 6935     Cidade de deus Distrito  District      <NA>      <NA>
## 6936 Complexo do Alemao Distrito  District      <NA>      <NA>
## 6937         Copacabana Distrito  District      <NA>      <NA>
## 6938          Guaratiba Distrito  District      <NA>      <NA>
## 6939 Ilha do Governador Distrito  District      <NA>      <NA>
## 6940            Inhaúma Distrito  District      <NA>      <NA>
## 6941              Iraja Distrito  District      <NA>      <NA>
## 6942        Jagarepagua Distrito  District      <NA>      <NA>
## 6943              Lagoa Distrito  District      <NA>      <NA>
## 6944          Madureira Distrito  District      <NA>      <NA>
## 6945               Mare Distrito  District      <NA>      <NA>
## 6946              Meier Distrito  District      <NA>      <NA>
## 6947             Pavuna Distrito  District      <NA>      <NA>
## 6948              Penha Distrito  District      <NA>      <NA>
## 6949          Portuaria Distrito  District      <NA>      <NA>
## 6950              Ramos Distrito  District      <NA>      <NA>
## 6951           Realengo Distrito  District      <NA>      <NA>
## 6952       Rio Comprido Distrito  District      <NA>      <NA>
## 6953            Rocinha Distrito  District      <NA>      <NA>
## 6954      São Cristóvão Distrito  District      <NA>      <NA>
## 6955         Santa Cruz Distrito  District      <NA>      <NA>
## 6956       Santa Teresa Distrito  District      <NA>      <NA>
## 6957             Tijuca Distrito  District      <NA>      <NA>
## 6958        Vila Isabel Distrito  District      <NA>      <NA>

Nós podemos plotar estes objetos através da função da plot:

# mapa do estado
plot(estado_rio)

# mapa da cidade
plot(cidade_rio)

No entanto, podemos também transformar estes shapefiles em um data.frame e fazermos as figuras no ggplot2, através da função fortify.

# criando um objeto do tipo data.frame para o pais
pais <- fortify(brasil)

## Regions defined for each Polygons

# observando o objeto
head(pais)

##        long       lat order  hole piece id group
## 1 -67.16084 -9.982599     1 FALSE     1  0   0.1
## 2 -67.15184 -9.978957     2 FALSE     1  0   0.1
## 3 -67.13828 -9.973469     3 FALSE     1  0   0.1
## 4 -67.11916 -9.967169     4 FALSE     1  0   0.1
## 5 -67.10312 -9.961880     5 FALSE     1  0   0.1
## 6 -67.08989 -9.900072     6 FALSE     1  0   0.1

# criando um objeto do tipo data.frame para o estado
estado <- fortify(estado_rio)

## Regions defined for each Polygons

# observado o objeto
head(estado)

##        long       lat order  hole piece   id  group
## 1 -44.11069 -23.18097     1 FALSE     1 6683 6683.1
## 2 -44.11097 -23.18097     2 FALSE     1 6683 6683.1
## 3 -44.11097 -23.18125     3 FALSE     1 6683 6683.1
## 4 -44.11153 -23.18125     4 FALSE     1 6683 6683.1
## 5 -44.11153 -23.18153     5 FALSE     1 6683 6683.1
## 6 -44.11180 -23.18153     6 FALSE     1 6683 6683.1

# criando um objeto do tipo data.frame para a cidade
cidade <- fortify(cidade_rio)

## Regions defined for each Polygons

# observando o objeto
head(cidade)

##        long       lat order  hole piece   id  group
## 1 -43.36049 -22.84699     1 FALSE     1 6929 6929.1
## 2 -43.37932 -22.85029     2 FALSE     1 6929 6929.1
## 3 -43.38162 -22.83656     3 FALSE     1 6929 6929.1
## 4 -43.39827 -22.84034     4 FALSE     1 6929 6929.1
## 5 -43.39690 -22.84560     5 FALSE     1 6929 6929.1
## 6 -43.40036 -22.84342     6 FALSE     1 6929 6929.1

# um mapa do estado do Rio de Janeiro
ggplot() + 
  geom_polygon(data = estado, aes(x = long, y = lat, group = group, fill = id), colour = "black", show.legend = FALSE) +
  coord_equal()

# um mapa para a cidade do Rio de Janeiro
ggplot() + 
  geom_polygon(data = cidade, aes(x = long, y = lat, group = group, fill = id), colour = "black", show.legend = FALSE) +
  coord_equal()

Esta transformação com a função fortify, no entanto, não é sempre necessária!

ggplot() + 
  geom_polygon(data = estado_rio, aes(x = long, y = lat, group = group, fill = id), colour = "black", show.legend = FALSE) +
  coord_equal()

## Regions defined for each Polygons

---

Exercício 4

Carregue este shapefile do Brasil, e:
* Crie um mapa do Estado do Rio de Janeiro, mas adicione cor somente à cidade do Rio de Janeiro;
* Crie um mapa da Cidade do Rio de Janeiro, mas adicione cor somente ao seu Bairro/Região (ou então escolha um bairro ou região à sua escolha caso você não more na cidade do Rio de Janeiro).

---

Usando informações da interação com a internet para criar mapas simples

Além de termos a opção de baixar arquivos dos tipos raster e shapefiles, nós também temos a opção de fazer o download de arquivos e mapas diretamente da internet, e projetarmos informações sobre eles. Muitas das funções que fazem isto estão no pacote ggmap.

• Argumento source:
  • google;
    
  • osm;
    
  • stamen;
    
  • cloudmade.
    
• Argumento maptype:
  • terrain;
    
  • toner;
    
  • watercolor;
    
  • satellite;
    
  • roadmap;
    
  • hybrid;
    
  • …

fundao <- get_map(location = c(lon = -43.2259, lat = -22.85769), source = "stamen", maptype = "watercolor", zoom = 14)
plot(fundao)

rio <- get_map(location = c(lon = -43.1729, lat = -22.90685), source = "stamen", maptype = "toner", zoom = 12)
plot(rio)

rio <- get_map(location = c(lon = -43.1729, lat = -22.90685), source = "google", maptype = "satellite", zoom = 14)
plot(rio)

Uma outra opção para conseguir mapas da internet está no pacote RgoogleMaps, através da função GetMap

rio_terreno <- GetMap(center = c(-22.92, -43.35), zoom = 11, maptype = "terrain", RETURNIMAGE = TRUE)
PlotOnStaticMap(rio_terreno)

rio_mobile <- GetMap(center = c(-22.92, -43.35), zoom = 11, maptype = "mobile", RETURNIMAGE = TRUE)
PlotOnStaticMap(rio_mobile)

rio_hybrid <- GetMap(center = c(-22.92, -43.35), zoom = 11, maptype = "hybrid", RETURNIMAGE = TRUE)
PlotOnStaticMap(rio_hybrid)